
//include standard libraries
#include <sstream>
#include <iostream> 
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <string.h> 
#include <fstream>
#include <vector>
#include <stddef.h>
#include <stdio.h>
#include <time.h>
#include <sys/time.h>
#include <algorithm>
#include <cassert>
#include <thread>

#include <typeinfo>
                        
#include <boost/math/distributions/gamma.hpp>
#include <boost/config.hpp>
#include <boost/generator_iterator.hpp>
#include <boost/random/linear_congruential.hpp>
#include <boost/random/uniform_int.hpp>
#include <boost/random/uniform_real.hpp>
#include <boost/random/gamma_distribution.hpp>
#include <boost/random/variate_generator.hpp>
#include <boost/regex.hpp>
#include <boost/math/distributions/normal.hpp>
#include <boost/random/mersenne_twister.hpp>
#include <boost/random.hpp>
#include <boost/random/normal_distribution.hpp>
#include "armadillo"

using namespace arma;
using namespace std;
using namespace boost;

//define some types because they are long and used a few times
typedef mt19937 base_generator_type;
typedef variate_generator<mt19937&, gamma_distribution<> > gamma_generator; 
typedef variate_generator<mt19937&, normal_distribution<> > normal_generator; 

vec Mahalanobis(mat x, rowvec center, mat cov);
//vec mahalanobis_arma( mat x ,  rowvec mu, mat sigma );
vec dmvnorm ( mat x,  mat mean,  mat sigma);


class icr_mcmc;
class icr_mcmc_state;
class icr_rands;
class rmvnorm;
class riwish;

class icr_mcmc_state{ //stores the current state of the MCMC
public:
  long iternum;//number of the iteration
  double alpha;//DP precision parameter
  rowvec m; //mean vector for G0
  mat V; //covariance for G0 
  mat V_inv; //invser of covariance for G0 
  mat S; //scale matrix for inverse-wishart part of G0
  rowvec p; // stick-breaking weights
  rowvec l; //configuration variable for assigning each data point to each 
            //component
  rowvec njstar;
  vector<vector<long>> l_rev;
  long nstar;
  mat mu_star; //one row vector per component
  cube sigma_star; //stores covariance matrix for each component
  double nu; //DOF for sigma_star riwish draws
};


class icr_mcmc{
public:
   icr_mcmc_state* state;
   long rowctr;
   icr_rands* rands; //pointer to rands object for drawing random numbers
   double logpn; //remainder of stickbreakings. Used in alpha updating
   long N; //number of components
   long ndata; //number of data points
   long niter; //number of mcmc iterations
   mat data; // matrix holding all the data
   long ncol; //number of data columns, ie number of dimensions
   //rowvec p;
   double L;
   mat H;
   rowvec h,r;
   string filename, outputname;
   ofstream fout;
   //priors for m
   rowvec a_m;
   mat B_m, B_m_inv;
   //priors for V
   double a_V;
   mat B_V;
   //priors for S
   double nu, a_S;
   mat B_S, B_S_inv;
   //priors for alpha
   double a_a, a_b;
  void getNextLineAndSplitIntoTokens(istream& str);
  void parse_csv(void);
  void find_nlines(void);
  void find_ncol(void);
  void blockfind(void);
  void stickbreaking_updates(void);
  void stickbreaking_first(void);
  void set_priors(void);
  void set_initial_values(void);
  void block_find(void);
  void update_indicators(void);
  void update_alpha(void);
  void update_params(void);
  void update_mustar_sigstar(void);
  void write_output(void);
  void run(void); //actually run the mcmc
  icr_mcmc(long,long,string,string);
  ~icr_mcmc();
};


class rmvnorm {
//class for drawing from a multivariate normal
public:
   icr_rands* rands; //pointer to rands object for drawing random numbers
   rowvec mu; //mean vector
   mat sigma; //covariance matrix
   mat chol_sigma; // cholesky decomposition of covariance matrix
   rowvec draw(void); //draws from the distribution
   rmvnorm(rowvec,mat, icr_rands*); //constructor
   ~rmvnorm();//destructor
};

class riwish{
//class for drawing from an inverse wishart distribution
public:
  double nu; //degrees of freedom
  //gamma_generator *gam_gens;
  vector<gamma_generator> gam_gens;
  mat S; //scale matrix
  mat chol_S; //cholesky factorization of S
  icr_rands* rands; //pointer to rands object
  mat draw(void); //draw from the distribution 
  mat draw_noninv(void); //draw from the distribution 
 // variate_generator<mt19937*, gamma_distribution<> >  makeOwnGamma(double);
  riwish(double, mat, icr_rands*);//constructor
  ~riwish();//destructor
};

class icr_rands {
//class for controlling all of your random number getting needs
public:
   time_t seed; //the stored seed
   normal_generator *normgen; //pointer to generator for normal deviates
   base_generator_type *using_seed; // the seed that is actually being used
                                    // passed around properly to all generators
   variate_generator
        <base_generator_type, uniform_real<> >  *uni; //uniform 0-1 gen
   double rand01(){return((*uni)()); } //returns uniform random deviate
   double randnorm(){ return((*normgen)()); } //return normal deviate
   gamma_generator rgamma_gen(double); //makes a gamma_generator
   void initialize_random_numbers(void); //initializes all random number stuff
   rowvec rdirichlet(rowvec); //returns a dirichlet distributed variable
   rmvnorm mvnorm_gen(rowvec, mat); //rerturns a multivariate normal generator
   riwish riwish_gen(double, mat); //returns an inv-Wishart generator
   long random_index(rowvec);
   icr_rands(); //constructor
   ~icr_rands(); //destructor
};

class icr_parser{
public:
  string infile, outfile, datafile;
  long nelem;
  vector<double> quantiles;
  vector<vector<double> > grid_params, dimens;
  vector<long> activedims, inactivedims;
  vector<long> set;
  vector<double>  getNextLineAndSplitIntoTokens(istream&);
  mat mu_star;
  cube sigma_star;
  rowvec njstar;
  double alpha;
  long niter, N, ncol;
  double ndata;
  vec output, output_tmp;
  mat evalpts;
  uvec permutation;
  //cube sigma_star_permuted;
  //mat mu_star_permuted;
  //rowvec njstar_permuted;
  double iternum;
  rowvec values;
  void read_csv(void);
  void run_parse(void);
  void make_grid(void);
  void test_read(void);
  icr_parser(string,string,string);
  ~icr_parser();
};


#include "icr_rands.cpp"
#include "icr_rmvnorm.cpp"
#include "icr_riwish.cpp"
#include "icr_mcmc.cpp"
#include "dmvnorm.cpp"
#include "icr_parser.cpp"

